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We apply the variational theory for fermions at finite temperature and high density, developed in 
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I. INTRODUCTION 



Ah initio models of dense nuclear matter at finite temperature are crucial to the un- 
derstanding of supernovae evolution , composition [3j and cooling JJ of protoneutron 
stars, gravitational wave emission spectrum from neutron star mergers [5|, and the analysis 
of heavy ion collision experiments [6|. Nuclear matter at zero temperature has been studied 
extensively using a variety of theoretical methods (see, e.g., [7j for a review). In contrast, 
;he corresponding many body problem at finite temperature has received little attention 
8fl. For example, most computer simulations of supernovae explosions use phenomenological 
equations of state, like the Lattimer-Swesty equation of state [9j or the equation of state due 



to Shen et al. 



10 



111- 



Three decades ago Friedman and Pandharipande (FP) carried out a seminal calculation 
of the equation of state of hot dense nuclear and neutron matter using a variational theory 



12J. Since then only a few other calculations have been carried out. The methods emp 



include Bloch-de Dominicis diagrammatic expansion 

I i 

self consistent Green's functions [lg], Dirac-Brueckner theory 



13|| , extended Brueckner theory 



Hartree-Fock theory 18 
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oyed 



id 



lowest order variational 20] methods. 



171 ]. relativistic Brueckner- 
jeiturbation theory with low momentum interactions []jj], and 



Most variational theories of dense quantum fluids originate from Jastrow's original sug- 
gestion that a reasonable approximation for the wave functions of an interacting system 
can be obtained by writing them as a product of the wave functions of a noninteracting 
system ($) and a product of pair correlation functions (Fy) {21], 

tfccJjFy*. (1) 

i<j 

The pair correlation functions are found by minimizing the ground state energy at zero 
temperature or the free energy at finite temperature. Different generalizations of this basic 
variational method have been applied to the many body problem with varying degree of 



success and sophistication 22]. In the approach used here the simple pair correlation func- 
tion in Eq. (PQ) is replaced by a pair correlation operator to include the effects of non-central 
correlations directly in the variational wave functions. The calculation of the energy expec- 
tation values are subsequently more difficult. Nevertheless a reasonably accurate calculation 
can be done using a combination of the Fermi hypernetted chain summation and the single 



operator chain summation approximations 



23| (and references therein). The combination of 



the Fermi hypernetted chain summation method for central correlations, the single operator 
chain summation for non-central correlations and some improvements introduced in later 
works (see next paragraph) is collectively known as the variational chain summation method 

The equation of state of dense nucleonic matter at zero temperature has been calculated 
using the variational chain summation method by FP, then by Wiringa, Fiks and Fabrocini 



25j and by Akmal, Pandharipande and Ravenhall (APR) [2J, |2fjfl. Each successive calcu 



lation was a significant improvement over the preceding one with respect to the variational 
method, and the models of the nucleon-nucleon interaction and three nucleon interaction 



used. Recently the three 
culated exactly (VCS/3) 



3ody cluster in the variational chain summation method was cal- 



271 ] unlike earlier calculations including APR where only the two 



body cluster terms were calculated exactly and all the higher order terms were calculated 
approximately (VCS/2). However VCS/3 has not yet been generalized to finite temperatures 
nor is a full zero temperature equation of state available at this moment. Thus in this work 
we have used the VCS/2 method and henceforth we will refer to VCS/2 simply as variational 
chain summation. 

We believe that variational chain summation supplemented by realistic nucleon-nucleon 
interactions and three nucleon interactions is at present one of the most reliable methods 
for studying dense many body systems. It compares reasonably well with the experimental 
bounds on the equation of state set by heavy ion collision experiments 



281 ] and benchmark 



quantum Monte Carlo calculations in neutron matter [2j|. It is, therefore, worthwhile to 
investigate its generalization to nonzero temperatures and arbitrary proton fraction. 

In extending variational chain summation to finite temperatures one is faced with a tech- 
nical challenge. Since wave functions used in variational chain summation are not mutually 
orthogonal, when constructing a variational theory at finite temperature one encounters the 



so called orthogonality corrections 



30 



3jj. The orthogonality corrections are not unique 



(they depend on the method of orthonormalization chosen) and more importantly their cal- 
culation requires the evaluation of the off-diagonal matrix elements of the Hamiltonian and 
unity. At present there exists no accurate method to calculate these off-diagonal matrix 
elements, especially for wave functions with operator dependent correlations. 
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In the only previous application of the variational method to finite temperatures, viz. in 
FP, the orthogonality corrections were simply ignored. Recently we showed that there exists 
a choice for the orthonormalization procedure such that the orthogonality corrections to the 



32fl (hereafter I). This way the orthogonality 



free energy vanish in the thermodynamic limit 
problem can be circumvented and at least for thermodynamic quantities the variational chain 
summation method can be extended to nonzero temperatures without having to worry about 
orthogonality corrections. 

In this paper we apply the formalism developed in I to symmetric nuclear matter and pure 
neutron matter. The present calculations can be regarded as an improvement on the finite 
temperature calculations due to FP and/or an attempt to generalize the zero temperature 
calculations of APR to finite temperature. This work is a contriubtion to the ongoing 
effort to understand the properties of nucleonic matter under extreme conditions of density, 
temperature and asymmetry starting from the bare interactions amongst nucleons which can 
be constrained by the experimental data on the nucleon-nucleon scattering in vacuum and 
the binding energies of light nuclei. 

The Hamiltonian we have used for this calculation is the same as that in APR. It consists 
of the Argonne v 18 nucleon-nucleon interaction 



interaction 



34 



331 ] . the Urbana IX model of three nucleon 



351 ] and leading order relativistic boost interactions [36]. For completeness 



we will discuss the Hamiltonian in some detail in the next section. In Section [111] we outline 
the variational chain summation method for finite temperature. First a short summary of 
the relevant parts of I is given. Thereafter the practicalities of the search for the variational 
minima are described. Section [IV] is devoted to our results. We first discuss the spin-isospin 
correlations in symmetric nuclear matter and pure neutron matter, and show that they are 



2J|. The fate of the neutral 



enhanced at finite temperature as they are at zero temperature 
pion condensation at finite temperature is discussed. The behavior of the nucleon effective 
mass, free energy, pressure and symmetry energy is also discussed in this section. Some of 
the limitations of our calculations are discussed in Section [V] We summarize our results in 
Section | 
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II. THE HAMILTONIAN 



The non-relativistic Hamiltonian used in the present calculations can be written in real 

space as 



» = E -t [(£ + i) + " r -] v - + E + H«) + £ W ( 2 ) 

i ij ijk 

where the kinetic energy operator takes into account the difference between the mass of a 
proton m p and the mass of a neutron m n . 

The Argonne V\% nucleon-nucleon interaction Uis,tj has the form 

18 
p=l 

The electromagnetic part v em is omitted from all nuclear matter studies including the present 
one. The strong interaction part has fourteen isoscalar operator terms 



O 



p=l,U 
ij 



1, cr, ■ o-j, Sij, (L • S) y , L?-, Ir?-<r< • (Tj, (L • S) 4 2 . ® [1, r< • r^] . 



(4) 



By convention the operators with even p < 14 have the Tj • Tj term while the ones with odd 
p do not. The three isotensor operators (p = 15, 16, 17) are given by 



O 



,p=15,18 



(3T zi T zj - Ti ■ Tj) ® (1, <Tj • (Tj, S, 



'.I i 



(5) 



And finally, the isovector operator (p = 18) is 



O 18 



(r zi + r zj ) . 



(6) 



The Argonne Vi 8 nucleon-nucleon interaction along with the CD BONN 



Nijmegen models 



37 



381 ] and the 



391 ] constitute the set of 'modern' phase-shift equivalent nucleon-nucleon 



interactions. These models fit the Nijmegen data base of proton-proton and neutron-proton 
scattering phase shifts up to 350 MeV with a x 2 /N data ~ 1. All of them include the long 
range one pion interaction potential but have different treatments of the intermediate and 
short range parts of the nucleon-nucleon interaction. 

The isotensor and isovector parts of Ui8,y, and the isovector part of the kinetic energy, 
are very weak and we will treat them as first order perturbations. In first order, these terms 



do not contribute to the energy of symmetric nuclear matter, which has total isospin T = 0. 
In pure neutron matter the isovector and the isotensor terms can be absorbed in the central, 
spin and the tensor parts of the nucleon-nucleon interaction. 

The Urbana IX model of three nucleon interaction V mXj ij k has two terms 

V vm>ij k = V?l + V? k . (7) 

The first term represents the Fujita-Miyazawa two-pion exchange interaction 

Vijk = Yscyc ^2tt ({T i ■ Tj,Ti ■ T k } {Xij, X ik } + \ [t \ ■ T j , T j • T k ] [Xij, X ik \) , (8) 

Xij = SijT n (rij) + <J, ■ a-jY^rij) , (9) 

with strength A 2n . The functions T n (rij) and K-( r ij) describe the radial shapes of the 
one-pion exchange tensor and Yukawa potentials. The term denoted by V^ k is purely phe- 
nomenological, and has the form 

Vgk = U y E7${r 1d )T$(r a ). (10) 

eye 

This term is meant to represent the modification of NA- and AA-contributions in the two- 
body interaction by other particles in the medium. The two parameters A 2n and Uq are 
chosen to yield the observed energy of 3 H and the equilibrium density of cold symmetric 
nuclear matter, p = 0.16 fm~ 3 . 

It is possible to incorporate the leading order effects of relativistic boost corrections as 



an interaction term into a non-relativistic Hamiltonian 36]. They give rise to two body 
interactions Sv^. As a first approximation, only the static part of the boost interaction is 
kept in the present calculation. Studies of light nuclei using the Variational Monte Carlo 



method 40] find that the contribution of the two-body boost interaction to the energy is 
repulsive, with a magnitude which is 37% of the V^ k contribution. However, the parameters of 
the Urbana IX model of three nucleon interaction were fixed without taking into account the 
boost interactions. In the new model of three nucleon interaction called UIX* the parameters 
are obtained by fitting the binding energies of 3 H and 4 He, and the equilibrium density of 
symmetric nuclear matter, includin g Sv jj. The strength of V^ k is 0.63 times that of V$ k in 
UIX, while VM remains unchanged 261 ] . 
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III. VARIATIONAL CHAIN SUMMATION AT FINITE TEMPERATURE 



A. The variational theory at finite temperature 

In variational calculations one assumes that a good approximation for the eigenstates of 



the interacting system of fermions is given by the correlated basis states 



12 



#i [7li(k, x)} = ^ g ( U ) 

K(k, x)\ (s (n Kj ^)) ®i K( k ' x )] 

where $j are the eigenstates of a free Fermi gas with occupation numbers (k, x) for single 
particle states with momentum k and x standing for any other quantum numbers of the 
single particle states (including spin and isospin). The occupation numbers rii (k, x) can 
take the values and 1. The pair correlation operator J 7 ^ encodes the effects of interactions. 
In the present calculations it has the following form 

8 

^ • (12) 

i=l 

Since the operators Ofj do not commute, the product of pair correlation operators need to be 
symmetrized with the symmetrization operator S to make the full wave function antisymmet- 
ric. Henceforth we will denote the static operator channels p = 1 — 6 interchangeably with 
the more descriptive symbols c, cr, a, or, t and tr which stand for central, central-isospin, 
spin, spin-isospin, tensor and tensor-isospin respectively. 

The ansatz Eq. fTTTT l is consistent with Landau's theory of Fermi liquids where it is assumed 
that the eigenstates of an interacting system have a one-to-one correspondence with the 
eigenstates of a non-interacting system. In this approach we assume that the mapping is 
accomplished by the correlation operator. In the spirit of Landau's theory we will hereafter 
refer to the rii(k,x) as the quasiparticle occupation numbers. 

An upper bound for the free energy of an interacting system can be obtained by using 



the Gibbs-Bogoliubov variational principle 



4l| 



F <F v = TY(p v H) + TTi(p v \np v ) , (13) 

where T is the temperature of the system and p v is a variational density matrix (not to be 
confused with the density p). The inequality is replaced by an equality when p v is the exact 
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density matrix of the system. Since we assume that the correlated basis states provide a 
good approximation to the eigenstates of the interacting system we want to construct p v 
from the correlated basis states. However, the correlated basis states, by construction, are 
not mutually orthonormal, i.e. in general 

^0 for i^j . (14) 



Hence, we need to orthonormalize them before they can be used 1 [3lfl. The orthonormal- 
ization procedure is not unique. Let |Oj) be the orthonormalized correlated basis states 
using one such orthonormalization method. Now, we can choose the following form for the 
variational density matrix, 

Pv = exp {-pH v )/Ti [exp {-(3H V ) ] (15) 
H v = E?\Qi)(ei\ , (16) 

where 3 = 1/T 

I — L 

and the variational spectrum can be approximated as a sum of quasiparticle energies [30|| 

El [n (k, x)\ = * (k, x; p, T) n (k, x) . (17) 

The quasiparticle energies e (k, x\ p, T) depend on the density p and the temperature T along 
with k and x. 

The entropy is given by 

S v (p,T) = -Tv(p v hip v ) 

= - In^XiPiT) ln {^( k ^;p> T )} 

+{1 - n(k, x- p, T)} ln{l - n(k, x; p, T)}} . (18) 

where n(k, x; p, T) is the mean quasiparticle occupation number of the single particle state 
(k, x) at density p and temperature T, 

n(k, x; p, T) = —— ^ , , \ ■ ( 19 ) 

exp [8{e (k, x; p,T) - p (p, T)}\ + 1 



1 It is possible to define the variational density matrix with the non-orthogonal correlated basis states. But 
in that case, the trace operation will involve their dual vectors and the free energy expectation value will 
still contain off-diagonal matrix elements. 
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Here jl is an effective chemical potential (see next subsection) and is fixed by 

A = ^n(k,a;;p,T) , (20) 

where A is the total number of particles in the system. We have set the Boltzmann constant 
ks to 1. 

The thermodynamic average of the Hamiltonian is given by 

Tv(p v H) = (H) OCB3jfi (k,x;p,T) (21) 

= (H) CBS ,n(k, X ;p,T) + E oc (p,T) (22) 

where (-H)cBs,fi(k ) x;/>,T) and (H) ocBs,n(k,x;p,T) are the expectation values of the Hamiltonian in 
the correlated basis states and orthonormalized correlated basis states, respectively, with the 
occupation numbers set to n (k, x; p,T). The orthogonality corrections are denoted by E oc . 

The calculation of the expectation values of various operators, especially the Hamiltonian, 
in the basis of correlated basis states is a nontrivial problem. The diagonal matrix elements 
and thus the expectation value in Eq. (l22l can be calculated by expanding it in powers of 
Tf, — 1. Schematically we have 



E v (p,T)/A = (H)aBa,nQve;p,T)/A 
2m 



h 2 

7j— k lv + dia g rams O> ^ ti) (23) 



where m is the average bare mass of a nucleon and k\ N is the mean square momentum per 
particle. The diagrams are many body integrals involving the potential v, the correlation 
operator T and the finite temperature Slater function It, 



Z r (r; p, T) = i ^ e lk r n(k, x; p, T) . (24) 

k,a; 

Within our scheme the pair correlation operator T and the Slater function It(?", p, T) encode 
all the relevant microscopic information about the many particle system at a temperature T 
and density p. 

Large classes of diagrams can be resummed using the variational chain summation 



methoc 
to SJ, 



For a review of the variational chain summation method the reader is referred 



24 



25 



26l | (and references therein). Here, we merely wish to point out that the 



variational chain summation method can be used to calculate the expectation values for any 
set of mean occupation numbers n (k, x). At zero temperature this set is a step function, 
but at finite temperatures states above the Fermi surface become populated. 

In contrast to E v (p, T) whose calculation involves only the diagonal matrix elements of the 
Hamiltonian in the correlated basis states basis, the orthogonality corrections E oc involve 
off-diagonal matrix elements of the Hamiltonian and the unit operator in the correlated 
basis states basis and they cannot be calculated within variational chain summation. In 
fact, there exists no method to calculate the off-diagonal matrix elements accurately for 
operator dependent interactions and correlation functions like the ones used in this study. 

In FP it was assumed that these orthogonality corrections are small. Recently it was 
shown in I that exploiting the fact that only a very small subset of states contributes to 
the trace in Eq. (1T31) , it is possible to construct an orthonormalization scheme where the 
orthogonality corrections to the free energy vanish in the thermodynamic limit 

-jp -> as A -> oo . (25) 

This proof does not rely on the detailed nature of the pair correlation functions, instead it 
follows from some general properties that any reasonable pair correlation function (including 
ours) is expected to possess. 

Hence using Eqs. (fl~8l 1221 1231 [25]) it is possible to rewrite Eq. f TT3T ) exactly as, 

F(p, T) < F v {p, T) = E v {p, T) - TS v (p, T) (26) 

In the above equation S v can be obtained trivially from Eq. (PT8l) once the single particle 
spectrum is known and E v is calculated using variational chain summation generalized to 
finite temperature as outlined above. For further details the reader is referred to I. 

B. The Optimization Procedure and the Pair Correlation Functions 

The variational free energy F v is optimized by varying both the single particle spectrum 
e(k, x; p, T) and the correlation operator Tij at all densities and temperatures. This should be 
contrasted with the calculations of FP where for densities p > 0.04 fm -3 the zero temperature 
correlation operator was used at all temperatures. We find that the temperature dependence 
of the correlation operator is weak but non-negligible (see later). 
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We parametrize the single particle spectrum e(k, x;p, T) by a simple effective mass ap- 
proximation 

^^-2^f)- (27 » 
In general it is possible for e(k, x; p, T) to have higher order terms in k. However, in our 
calculations F v was found to be insensitive to any such dependence. It is also possible for 
e(k, x; p, T) to contain a momentum independent term which depends on p and T only, 



u(p, T), but such a term will be absorbed in the definition of p [42 1 



p(p,T) = p(p,T)-u(p,T) (28) 

where p is the true chemical potential of the system. 

As mentioned earlier, in variational chain summation the expectation values of various 
operators are calculated by expanding the matrix elements in powers of Tfj — 1 and then 
resumming large classes of terms. The correlation operator JF^ is itself calculated by solving 
the Euler-Lagrange equations which are obtained by minimizing the lowest order terms ( 
viz., the sum of the kinetic energy term and the two body cluster term) in this expansion 
for the Hamiltonian, with the nucleon-nucleon interaction Vij replaced by Vij — Xij, where 



Vij 



J2 ^(ryJO?-, (29) 

p=l,14 

A« = E XP ^)°%- ( 3 °) 

p=l,8 

The variational parameters a p are meant to simulate the quenching of the spin-isospin in- 
teraction between particles % and j, due to flipping of the spin and/or isospin of particle i or 
j via interaction with other particles in matter. We use 

c? = 1 for p = 1 and 9 (31) 
a p = a otherwise . (32) 

The correlation functions f p {r) are made to satisfy the additional healing conditions 

f p (r > (F) = S pl (33) 
df p 

= (34) 

dP 



dr 
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where 



d p = d t for p = 5, 6 
d p = d c for p ^ 5, 6. 



(35) 
(36) 



The above constraints completely determine the A p (r). 

The Euler-Lagrange equations are solved in the spin S and isospin T channels. As an 
example, in the S = channel the Euler-Lagrange equations become 

h 2 

- ^ [0r,5=oV 2 /r,s=o + 2V0t,s=o • V/r,s=o] + (v T ,s=o - A r ,s=o) h,s=o<pT,s=o = (37) 
where 



[l-(-l)< T+s >#(r;p,T 



1/2 



(38) 



The relationship between the potential and the correlation functions in the T, S channels 



and those in the operator channels p = 1 — 8 are given in [23]. The Euler-Lagrange equa- 
tions in the S = 1 channels are considerably more complicated because the contribution of 
the tensor and the spin-orbit terms give rise to three coupled differential equations for the 



correlation functions in each isospin channel. For more details the reader is referred to 23]. 

To illustrate the relative effect of density and temperature on the pair correlation functions 
we show in Figs. Q] and [2] the pair correlation functions in the central (f c ) and the tensor- 
isospin (f tT ) channels for p = p and 1.5p and T = and 10 MeV. As we mentioned in 
the beginning of this section, thermal effects (at least in the range of temperatures we are 
interested in) are not negligible but they do not change the behavior of the / p 's qualitatively. 

The variational free energy F v (p,T) thus becomes a function of the four variational pa- 
rameters to*, d c , d t and a. The optimal values for these parameters and hence the best value 
F v at each density and temperature is found by minimizing the constrained free energy 
defined as 



F V + AA 



(Ic - I) 2 + (gir + I) 2 



where 



I c = pj d 3 r(l-g c (r)) 



(39) 

(40) 
(41) 



i,j=l,N 
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and g c {r) is the pair distribution function. Laws of conservation of mass and charge demand 
that, 

Ic = 1 (42) 
I T = -3 . (43) 

Of course, for pure neutron matter only Eq. (1421 ) is applicable and hence, only the first term 
within square brackets in Eq. ( 1391) is kept. 

In the right hand side of Eq. f[39l) penalty term is added to F v to make sure that these laws 
of conservation of mass and charge are approximately satisfied in the variational calculations. 



In our calculations we choose A to be 1000 MeV [2J|. This keeps I c and I T to within about 
5% of their exact values for all values of p and T. For most values of p and T the conservation 
laws are satisfied at the level of one percent or less. 

The temperature dependence of the correlation operator is important in satisfying the 
sum rules Eqs. (|42~I) and f|43l) to a reasonable degree of accuracy. For example, at p = p and 
T = 10 MeV, the deviation of I c and I T in our calculations from their true values (Eqs. (1421) 
and (1431 ) is ~ 0.04% and ~ 0.3%, respectively. For comparison, employ the methodology 
used by FP, the so called 'frozen correlation' method, where the correlation operator is 
taken from the zero temperature calculations, only m* is varied at finite temperature, and 
the minimization is carried out with F v only. At the same values of p and T the deviation 
of I c and I T from their correct values is ~ 13% and ~ 20%, respectively, in this case. 

In practice we varied m*/m, d t /r 0l d c /r and a during the search. Here r , is the unit 
radius, defined such that 

4 

-irrlp = 1 . (44) 

Due to technical reasons involving the spacing on the grid on which the integrations are done 
in our computer program, d c /r is varied in fixed steps while the other three parameters are 
allowed to vary continuously. In particular, 

Ic = kj t (45) 

could only take integer values where 2l t is the size of the grid on which most of the integrations 
are carried out. In our calculations we set l t to be 64 while this was 32 in the calculations due 
to APR. We also made improvements in the way the integrations are done while calculating 

13 



the chain functions. While in APR a simple midpoint Euler method was used, we employed 
Gaussian quadrature for the same purpose. 

The actual variational search for each value of L is carried out using a downhill simplex 

n 

routine [43fl. Within the context of variational chain summation the simplex search algorithm 



25l |. Finally, the best parameters 



for the search parameters was first implemented in Ref. 
and F con are found from a quadratic fit to the corresponding values at the three best values of 
l c . This is an improvement over APR where a simple grid search was done in the parameter 
space. 

Typically the changes in energy due to these improved numerics are small (~ 0.5 MeV 
per particle at p = p and T = 0), but the short range parts of the chain functions become 
much more well behaved as a result of these improvements, i.e., the quality of the final wave 
functions obtained is much improved 2 . 



IV. RESULTS 



A. Two Body Densities and Spin-isospin Correlations 

The effect of nuclear interactions on two particle correlations in medium can be inferred 
from the two body densities p^(r), defined such that, 

( £ B( riJ )0^ T , p = A I d\B(r)p p 2 (r) , (46) 

where (• ■ -)t, p denotes the thermal average at temperature T and density p. We show the 
static two body densities in Fig. [3] for T = 0, 10 and 20 MeV at saturation density p = p in 
symmetric nuclear matter. Unsurprisingly, they have the expected asymptotic behavior 

p\{r -> oo) = p5 lp . (47) 

The central two body density p\ is particularly interesting because it is p times the probability 
of finding two particles separated by a distance r. In a non-interacting gas all non-central 
(p > 1) two body densities vanish; their large magnitude in nuclear matter is due to strong 
spin-isospin correlations introduced by the nuclear interaction. 

2 Some of the corresponding improvements for the zero temperature version of the computer program were 
first carried out by Jaime Morales 
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The effect of increasing temperature on the two body densities, within the context of our 
calculations, is twofold. Firstly, increasing temperature weakens the effects of Pauli blocking 
in the Euler-Lagrange equations (like Eq. (1371) ) for the pair correlation functions . This, tends 
to enhance short range correlations among nucleons. On the other hand, when a nonzero 
temperature is introduced the integral equations used in variational chain summation to 
calculate the the two body densities from the pair correlation functions, it tends to suppress 
the spin-isospin correlations. In the regime of temperatures and densities discussed in this 
work, the magnitudes of the two effects are comparable and the net temperature dependence 
of the two body densities is a combination of the two. 

Generally particles are correlated over the longest range due to tensor interactions. In 
our calculations d t is one measure of this range. In Figs. 0] and [5] we show d t for symmetric 
nuclear matter and pure neutron matter, respectively, for various temperatures. In Fig. [4] 
we see a dramatic decrease in the value of d t /r below p ~ 0.06 fm~ 3 for T = 4 MeV in 
symmetric nuclear matter. This is a precursor to cluster formation. We will discuss this 
topic briefly in Section [V] 

One of the most interesting features of these plots is the sharp change in magnitude of d t /r 
at p ~ 2p for symmetric nuclear matter and at p ~ po for pure neutron matter. A similar 
feature was obtained in the zero temperature calculations of APR. In APR it was argued 
;hat this feature resulted from a first order phase transition due to neutral pion condensation 
Their argument was based on the fact that the energy per particle as a function 
of density showed a kink indicating a first order phase transition, and that the isovector 
spin longitudinal reponse showed an enhancement and softening in the high density phase 
compared to the low density phase indicating enhanced pion exchange interactions between 
nucleons in the high density phase. Recently there has been some indication of experimental 
evidence supporting the enhancement of pionic modes in nuclei (see [471 ] for a review). 

The isovector spin longitudinal static structure function is defined as 



1 



((OD 



\(Ol)tJ 2 ) 



(48) 



A 



where 



.4 




(49) 



i=l 



15 



Here rj are the positions of the nucleons and t is a unit vector in the isospin space. This quan- 
tity is also the sum of the isovector spin longitudinal dynamic response function i?z(q, a;), 

S , i(?) = 3/o°°i2L(q,w)da; (50) 
Rl&u) = J2i,jPi\(J\0 L (q)\I)\ 2 S(u;j - u: - u) (51) 

where I and J are eigenstates of the Hamiltonian with eigenvalues uj and uj, respectively, 
and pi is the probability of the state I occurring in the thermal ensemble. 

It is also possible to define a mean energy E L (q) of the isovector spin longitudinal response 

by, 

W L (q) = jj uoR L (q } uj)duj (52) 

s / \ W L (q) 

= ^ (53) 

Both Si(q) and Wi{q) and hence E L (q) can be calculated from the two body densities p^{r). 

In Fig. [6] we plot Sl and El in symmetric nuclear matter at T = 4, 10 and 16 MeV for 
various densities. For all three temperatures we see a big enhancement of Si at q ~ 1.3 fm -1 
due to stronger spin-isospin correlations at higher densities. However, for T = 4 MeV the 
enhancement develops quite suddenly at p = 2p while at T = 10 MeV and 16 MeV there 
is a smoother evolution of the enhancement as the density is increased. Similarly the mean 
energy E L develops a dip at q ~ 1.3 fm -1 for all three temperatures, but quite sharply at 
T = 4 MeV and relatively smoothly at 10 and 16 MeV as functions of density. 

The same quantities for pure neutron matter at T = 14, 22 and 30 MeV are plotted in 
Fig. d Again, we see a large enhancement in the value Si for p > p at q ~ 1.3 fm" 1 at 
all three temperatures. However the enhancement develops suddenly in the case of T = 14 
MeV, and relatively smoothly for T = 22 and 30 MeV. Also, the mean energy E L shows a 
dip for densities p > p at all temperatures, but this effect develops more smoothly at higher 
temperatures. 

Our calculations thus show that the isovector spin longitudinal response for symmetric 
nuclear matter is enhanced and softened at densities p > 2p at finite temperatures. But 
whereas for T < 10 MeV this enhancement and softening develops quite suddenly, there is 
a smoother evolution for T > 10 MeV. For pure neutron matter there is enhancement and 
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softening of the isovector spin longitudinal response at densities p > po, but it develops quite 
suddenly for T < 22 MeV and quite smoothly for T > 22 MeV. 

The difference between the expectation values of the pion number operator in a system 
of A interacting nucleons and in a system of A isolated nucleons is called the pion excess. 
The part (5^ (q)) of the pion excess at a momentum q, exclusively due to one pion exchange 



interactions, can be calculated from Si{q) 2d, 148] . In Figs. M and M we show (S%(q)) 
for symmetric nuclear matter and pure neutron matter. As with S^q), (8-p(q)) shows an 
enhancement at higher densities. For symmetric nuclear matter this enhancement develops 
suddenly for T < 10 MeV and smoothly for T > 10 MeV. In pure neutron matter T « 22 
MeV marks the temperature below which the enhancement is sudden and above which the 
enhancement is smooth. 

The above observations along with the results for the equation of state (described in 
Section HV CI below) leads us to conclude that the first order phase transition due to neutral 
pion condensation has a critical temperature of T c m 10 MeV for symmetric nuclear matter 
and T c 22 MeV for pure neutron matter. Above the critical temperature the enhancement 
in the spin-isospin correlations still persists at higher densities, but now there is a smooth 
crossover between the low density phase and the high density phase. 



B. The Nucleon Effective mass 



The effective mass m* is not only a variational parameter in our calculations, but also 
a quantity of considerable physical importance in determining the thermal properties of 
nuclear matter. Our results for m* as a function of density in symmetric nuclear matter 
and pure neutron matter for various temperatures is shown in Figs. [10] and [TH We wish to 
point out that the variational minimization in pure neutron matter is not very sensitive to 
variations of m*, especially at low temperatures and the high density phase. In particular 
the uncertainty in our estimate of rrf for the high density phase in pure neutron matter 
for the three lowest temperatures reported (4, 8 and 12 MeV) is comparable to the total 
temperature dependence of m* at these temperatures. 

The curves of m* vs p show sharp changes at p ~ 2po in symmetric nuclear matter and 
at p ~ po in pure neutron matter. The origin of these changes is the enhancement of the 
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spin-isospin correlations discussed earlier. 

Our calculations also show an enhancement of m* at low temperatures in symmetric 
nuclear matter and to a much lesser extent in pure neutron matter. This is clearer from 
Fig. [12] where m* has been plotted as a function of temperature for p = 0.5p , Po an d l-5p 
in symmetric nuclear matter. For example, for p ~ po, m* at T = 2 MeV is 0.85m and at 
T = 4 MeV is 0.81m while at T = 20 MeV it is 0.69m which is very close to the value 0.7m 
obtained from optical potential models of nucleon-nucleus scattering 

HQ- 
It is important to realize that the effective mass, defined as 

m* W = k«L , (54) 

can in general depend on the momentum k. The momentum independent effective mass 
m*(p, T) that is used in the variational calculations is an weighted average of this momentum 
dependent m*(k). It is difficult to establish the actual weighting given to each single particle 
state. However it is not difficult to convince oneself that at any given temperature and 
density the maximum contribution comes from states with single particle energies e(k) such 
that 

\e(k)-e(k F )\ <ttT (55) 

where kp is the Fermi momentum. 

At very low temperatures only those single particle states which are very close to the 
Fermi surface can contribute. Theoretical calculations and experimental evidence suggests 
that the (momentum dependent) m*{k) in symmetric nuclear matter has a big enhancement 

In our calculations enhancement in m* at 



at and near the Fermi surface 



50, 



51 



52 



low temperatures is probably due to this effect. 

In FP no such enhancement was seen. For comparison, in Fig. [12] we have also plotted 
the m* at p = po obtained by simply minimizing F v and using correlation operator from the 
T = calculations, i.e. the frozen correlation method. This method is closest to the one used 
in FP albeit with different nucleon-nucleon interaction and three nucleon interaction. And 
indeed, we do not see any significant variation of m* with temperature. In fact, the values of 
m*/m so obtained are 0.66, 0.68, 0.71, 0.73 at T = 5, 10, 15 and 20 MeV, respectively. These 
are very close to the corresponding values 0.65, 0.67, 0.70 and 0.73 obtained by FP where 
the Urbana v 14 model of the nucleon-nucleon interaction was used along with a density 
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dependent term to simulate the effects of the three body forces. 

It appears that we are able to capture some of the subtle correlations which influence 
the low energy quasiparticle excitations better because of the improvements introduced in 
our calculations. These improvements include the fact that in the variational search all the 
parameters were varied and the pair correlation functions were calculated at all temperatures, 
and the fact that we insist on satisfying the laws of conservation of mass and charge at the 
level of a few percents. However, this is simply one possibility since in general the free energy 
is not very sensitive to the variations in the effective mass and it is difficult to delineate the 
effect that each individual component in our method has on the effective mass. This topic 



deserves further investigation, and we hope to address it in more detail in future work 



55fl. 



C. The Equation of State 

As we discussed earlier the variational minimization for the optimal parameters d t , d c , a 
and m* is carried out with Vis + V mx only. Following APR some additional corrections were 
added to the energy. 

1. The relativistic boost corrections 5vb + (V* 1X — V mx ) which were discussed earlier. 

2. An estimate for the second order perturbative corrections AE 2 is given by bE 2 B the 
difference between the two body cluster energy obtained by minimizing the contri- 
butions in each partial wave and the two body cluster energy obtained from the Tij 
defined above. 

3. An additional ad-hoc correction term — 72p 2 e _73P , with 72 = 1996 MeV fm 6 and 73 = 
15.24 fm 3 is added to the symmetric nuclear matter energy to get the correct saturation 
energy at zero temperature. These values are slightly different from those in APR (72 = 
2822 MeV fm 6 and 73 = 18.34 fm 3 ) because our variational energies and variational 
parameters at zero temperature are slightly different from those in APR owing to the 
improvements in the method described earlier. 

Of the above only the relativistic boost corrections have a temperature dependence. The 
temperature dependence of the perturbative corrections was found to be negligible and was 
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neglected. Thus our best estimate for the free energy is given by 



tot 



F v + E c 



E v — TS V + E c 



(56) 
(57) 



where E c is the sum of the correction terms mentioned above. 

The free energy F tot calculated using the method outlined in the last section is shown 
in Figs. [13] and [EH for symmetric nuclear matter and pure neutron matter, respectively, 
for temperatures T = 4 — 30 MeV. The most striking feature in these figures is the sharp 
change in the slope of the free energy at low temperatures for both symmetric nuclear matter 
(T < 12 MeV, p ~ 0.30 fm~ 3 ) and pure neutron matter (T < 24 MeV and p ~ 0.18 fm~ 3 ) 
(see inset in Figs. [T3l and [T4l) . 

During the variational search the global minimum of the free energy as a function of the 
variational parameters jumps from one local minimum to another. In Fig. [Tj3 we show the 
density at which this transition from the low density phase to high density phase occurs in 
our calculations. This results in a first order phase transition. In light of our discussion 
in Section IIV Al we identify this phase transition with the phenomenon of the neutral pion 
condensation. As we mentioned in the earlier the critical temperature for this transition is 
T c « 10 MeV for symmetric nuclear matter and T c « 22 MeV for pure neutron matter. 

At subsaturation densities and low temperatures symmetric nuclear matter undergoes 
another first order phase transition, the liquid-gas transition. In our calculations the critical 
temperature for the liquid-gas phase transition is T c 21 MeV and the critical density is 
about p c ~ 0.3p - If the calculations are done using the variational free energy F v only, then 
the critical temperature for the liquid gas phase transition comes out to be about 16 MeV. 
In the calculation due to FP the critical temperature is about 17.5 MeV. 

The pressure can be calculated from, 



Our estimates for the free energy F(p,T) (F tot of Eq. (f57l)) and the pressure P(p,T) of 
symmetric nuclear matter at T = 10 MeV from the variational chain summation calculations 
have been plotted in Figs. [16] and [17] along with the same from some of the other popular 
equations of state viz. the equation of state due to FP, the Lattimer-Swesty (LS) liquid 



P(P, T) = p 



]2 OF I A 
dp 



T 



(58) 
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droplet model calculations with incompressibility K = 180 and 220 MeV [9] and the equation 
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of state due to Shen et al. (STOS) [10, Qj| which is based on relativistic mean field theory. 



The inset shows the pressure in the low density region in greater detail. We have kept the 
thermodynamically unstable region in the P vs p curve to facilitate comparison with the 
other equations of state. Due to the enhanced spin-isospin correlations our equation of state 
is relatively soft at p ~ 2p fm~ 3 but hardens at higher densities. 

Various astrophysical phenomena depend sensitively on the symmetry energy E sym [lj and 



the symmetry free energy F S} 



561 ] which are defined as 



1 d 2 E 



2 dS 2 
1 d 2 F 



2 d5 2 



8=0 



8=0 



Here 5 is the asymmetry parameter given by 



8 



Pn - Pp 



(59) 
(60) 

(61) 



where p n is the neutron density and p p is the proton density. To the leading order in S, E syin 
is given by the difference between the energy in pure neutron matter E PNM and the energy 
in symmetric nuclear matter E SNM . Similarly F sym is given by the difference between the free 
energy in pure neutron matter -F PNM and symmetric nuclear matter -F SNM - 



E F 



E s 



■ + 0(5 2 ) 
+ 0(5 2 ) 



(62) 
(63) 



In Figs. [El and [19] we show E PNM — E SNM and F PNM — F SNM as functions of density at various 
temperatures. 



V. DISCUSSION 



In this section we will discuss the limitations of our calculations. First, let us consider 
the validity of the quasiparticle picture at finite temperature. This assumption is motivated 
by Landau's theory of Fermi liquids. However Landau's theory in its original formulation 
applies only to systems at zero or very low temperatures. It is justified on the basis of the 
fact that low lying quasiparticle excitations (as defined by the poles in the single particle 
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propagator) are long lived at low temperatures [57|, [58||. At higher temperatures this is not 
true and one cannot associate well defined quasiparticle excitations to the system. 

It is, however, important to recognize that there are two separate issues involved here. 
On the one hand, we have the thermodynamical statement that the free energy of the system 
can be written as a functional of quasiparticle occupation numbers and that in particular 
the entropy is given by Eq. (fT8l ). On the other hand, we have one particular microscopic 
justification of the previous statement in terms of the long lived nature of the quasiparticle 
excitations. It has been shown in the past that although the latter is true only for very low 
iemperatures, the former, under very general conditions, is true for arbitrary temperatures 



59 



60 



61 



62], Q- With this in mind we think that our final estimate for the free energy 



is more reliable than the individual components that entered the calculations, e.g. the 
assumption of an one-to-one mapping between the states of a non interacting system and 
an interacting one is not exactly true at for the states relevant at any finite temperature. 
However, even with this caveat in mind, we believe that our predictions about the microscopic 
structure of the system are useful in accounting for at least the gross features if not the details 
in the many body system. 

Our description of the high density phase is clearly incomplete. Our wave function does 
not include any true long range order. It has been shown from very general arguments 
by Landau and Peierls that long range order where the order parameter varies only in one 



dimension is unstable at finite temperatures 



64 ]. However, condensates that vary in two 



or three dimensions can exist at any temperature. If the high density phase is indeed a 
state with broken symmetry then it is possible that the first order transition changes into a 
second order transition beyond critical temperature (instead of a smooth crossover from the 
low density phase to the high density phase seen in our calculations). From our calculations 
we cannot make any prediction about this possibility. 

In the low density and low temperature region for symmetric nuclear matter the uniform 
fluid is no longer the optimal configuration and few body clusters are expected to emerge 
(see e.g. ]65]). We see the precursor to deuteron clustering in our calculations, e.g. in Fig. EO 
we see that the pair distribution function for symmetric nuclear matter at p = 0.02 fm -3 
and T = 2 MeV shows a pronounced bump at r ~ 1 fm. Genuine clustering of three or more 
particles, however, is not possible in our variational chain summation calculations because 
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our wave functions are too simple to describe such phenomena. 

Another phenomenon which occurs at very low densities and temperatures in nuclear 



matter is superfluidity due to the formation of Cooper pairs 



66 



67|] . Pairing, especially 



in very low density pure neutron matter and neutron star matter, has been the subject of 



extensive research in the recent past 



681 ] (and references therein). The wave functions that 



we have used do not include the effects of pairing in the long range part. However we would 
like to point out that it is possible to develop a variational theory based on correlated basis 
states for superfluid systems at zero temperatures It is quite possible that this 

theory can be extended to finite temperatures in a manner very similar to the one used here 
for normal systems. 

At very high temperatures (T > 50 MeV for symmetric nuclear matter) thermal excita- 
tions of the pionic degrees of freedom, the 'thermal pions' [zi] become important and should 
be explicitly included in the equation of state {3] . Our present calculations do not have 
any explicit non-nucleonic degrees of freedom, thus they cannot be extended to very high 
temperatures without the explicit inclusion of the thermal pions. The calculation of the 
pionic spectrum in nuclear matter is a challenging task. However, as the first step, the sim- 
plified quasiparticle treatment of the pions along with the full variational chain summation 



calculation for the nucleons as done in Ref. 



72J might prove useful. 



Finally, we would briefly discuss the case of arbitrary proton fraction. Pure variational 
chain summation can only be used to calculate the equation of state of symmetric nuclear 
matter and pure neutron matter, and not for asymmetric nuclear matter. This problem can 
be circumvented by using the so called quadratic approximation, where one assumes that 
the interaction energy at an arbitrary asymmetry 5 has a purely quadratic dependence on 
S. Using this fact the energy for any S can be found from the results at 6 = (symmetric 



26j, 
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nuclear matter ) and S = 1 (pure neutron matter) 

It has been shown in Ref. [73j that within the scheme of variational calculations (although 
with a lower order approximation for the non-central correlations as compared to the varia- 
tional chain summation method) the quadratic approximation is quite accurate, at least up 
to the leading order terms in the interaction energy. This calculation was done for the zero 
temperature low density phase, with the Urbana v 14 model for nucleon-nucleon interaction 
and a density dependent model for three nucleon interaction. Thus far, the calculation has 
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not been repeated for finite temperatures, the high density phase or for the modern sets of 
nucleon-nucleon interaction and three nucleon interaction. It is, however, plausible that the 
result is true for the aforementioned conditions as well. 



VI. CONCLUSION 

The calculation of the properties of nuclear matter in supernovae, neutron stars and 
heavy ion collisions starting from the properties of nucleons in vacuum and light nuclei is an 
outstanding problem in nuclear many body theory. In a preceding paper (I) we showed that 
when the variational theory is extended to finite temperatures, the free energy and the single 
particle energies do not have any corrections due to the non-orthogonality of the correlated 
basis states used in these variational theories. In this paper we showed how this formalism 
can be used for practical calculations in hot dense nucleon matter. 

We calculated the equation of state of symmetric nuclear matter and pure neutron matter 
and discussed the fate of the neutral pion condensation at finite temperature. The first order 
phase transition was seen to have a critical temperature of about 10 MeV for symmetric 
nuclear matter and about 22 MeV for pure neutron matter. However, the the enhancement 
of the spin-isospin correlations and as a result the softening of the equation of state, near 
p ~ 2p for symmetric nuclear matter and near p ~ Po f° r pure neutron matter, persists 
even beyond the critical temperature. We also discussed the behavior of the nucleon effective 
mass and saw an enhancement at lower temperatures. We proposed one possible explanation 
for this effect. 

The pair correlation functions Tij generated during our variational minimizations can be 
used to calculate the equation of state of asymmetric nuclear matter within the quadratic 



approximation 



interactions 



76 



the single particle energies in dense nuclear matter 0, 



12 

our future work |55 | 



2f|, effective 



and transport properties [7_S|]. Some of these topics will be explored in 
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Figure 1: (Color online) The pair correlation function in the central channel of symmetric nuclear 
matter various densities and temperatures. 
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Figure 2: (Color online) The pair correlation function in the tensor-isospin channel of symmetric 
nuclear matter various densities and temperatures. 
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Figure 3: (Color online) The two body densities in symmetric nuclear matter at saturation density 
and different temperatures. 
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Figure 4: (Color online) The tensor correlation length in symmetric nuclear matter. The numbers 
alongside the curves denote the temperature in MeV. 
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Figure 5: (Color online) The tensor correlation length in pure neutron matter. The numbers 
alongside the curves denote the temperature in MeV. 
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Figure 6: (Color online) The isovector spin longitudinal static structure function and mean energy 
in symmetric nuclear matter. 
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Figure 7: (Color online) The isovector spin longitudinal static structure function and mean energy 
in pure neutron matter. 
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Figure 8: (Color online) The pion excess due to one pion exchange interactions in symmetric nuclear 
matter. 
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Figure 9: (Color online) The pion excess due to one pion exchange interactions in pure neutron 
matter. 
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Figure 10: (Color online) The effective mass in symmetric nuclear matter. The numbers alongside 
the curves denote the temperature in MeV. 
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Figure 11: (Color online) The effective mass in pure neutron matter. The numbers alongside the 
curves denote the temperature in MeV. 
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Figure 12: (Color online) The effective mass in symmetric nuclear matter at p = 0.5po, Po and 
1.5po- Also plotted is the effective mass at p = po, obtained by using the frozen correlation method 
described in the text. 
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Figure 13: (Color online) The free energy in symmetric nuclear matter. The numbers alongside 
the curves denote the temperature in MeV. The inset shows in more detail the region of transition 
from the low density phase to the high density phase. 
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Figure 14: (Color online) The free energy in pure neutron matter. The numbers alongside the 
curves denote the temperature in MeV. The inset shows in more detail the region of transition from 
the low density phase to the high density phase. 
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Figure 15: The transition density for neutral pion condensation for symmetric nuclear matter and 
pure neutron matter in the variational calculations. 
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Figure 16: (Color online) A comparison of the free energy obtained using various methods at T = 10 
MeV. For an explanation of the various curves see the text. 
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Figure 17: (Color online) A comparison of the pressure obtained using various methods at T = 10 
MeV. The inset shows the low density region in more detail. For an explanation of the various 
curves see the text. 
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Figure 19: (Color online) A comparison of F PNM — F SNM at various temperatures 
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Figure 20: (Color online) The normalized pair distribution function at T = 4 MeV and p = 0.02 
fm -3 shows the effects of deuteron clustering. The same function has been plotted for T = 20 and 
p = 0.02 fm -3 and T = 4 MeV, and p = 0.16 fm -3 for comparison. 
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